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ABSTRACT 

We present a comprehensive synthesis model for the AGN evolution and the growth of super- 
massive black hole in the Universe. We assume that black holes accrete in just three distinct 
physical states, or "modes": at low Eddington ratio, only a radiatively inefficient, kinetically 
dominated mode is allowed (LK); at high Eddington ratio, instead, AGN may display both a 
purely radiative (radio quiet, HR), and a kinetic (radio loud, HK) mode. We solve the continu- 
ity equation for the black hole mass function using the locally determined one as a boundary 
condition, and the hard X-ray luminosity function as tracer of the AGN growth rate distri- 
bution, supplemented with a luminosity-dependent bolometric correction and an absorbing 
column distribution. Differently from most previous semi-analytic and numerical models for 
black hole growth, we do not assume any specific distribution of Eddington ratios, rather we 
determine it empirically by coupling the mass and luminosity functions and the set of fun- 
damental relations between observables in the three accretion modes. SMBH always show 
a very broad accretion rate distribution, and we discuss the profound consequences of this 
fact for our understanding of observed AGN fractions in galaxies, as well as for the empirical 
determination of SMBH mass functions with large surveys. We confirm previous results and 
clearly demonstrate that, at least for z < 1.5, SMBH mass function evolves anti-hierarchically, 
i.e. the most massive holes grew earlier and faster than less massive ones. For the first time, 
we find hints of a reversal of such a downsizing behaviour at redshifts above the peak of 
the black hole accretion rate density (z w 2). We also derive tight constraints on the (mass 
weighted) average radiative efficiency of AGN: under the simplifying assumption that the 
mass density of both high redshift (z ^ 5) and "wandering" black holes ejected from galac- 
tic nuclei after merger events are negligible compared to the local mass density, we find that 
0.065 < Co(eiad} < 0.07, where is the local SMBH mass density in units of 4.3 x IO^Mq 
Mpc^'^. We trace the cosmological evolution of the kinetic luminosity function of AGN, and 
find that the overall efficiency of SMBH in converting accreted rest mass energy into kinetic 
power, Ekin, ranges between ekin — 3 ~ 5 x lO^'^, depending on the choice of the radio core 
luminosity function. Such a "kinetic efficiency" varies however strongly with SMBH mass 
and redshift, being maximal for very massive holes at late times, as required for the AGN 
feedback by many galaxy formation models in Cosmological contexts. 
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1 INTRODUCTION 

Soon after their discovery, it was realized that Quasars (QSOs) were 
a strongly evolving Cosmological population, that could effectively 
be used as a tracer to study the structural properties of the evolving 
Universe ( [Longair 1966) . The physical understanding of QSOs and 
of their lower-luminosity counterparts (both generally called Ac- 
tive Galactic Nuclei, AGN) as accreting supermassive black holes 
(SMBH) immediately led to speculations about the presence of 
their dormant relics in the nuclei of nearby galaxies, with Soltan 
(1982) first proposing a method to estimate the mass budget of 



SMBH based on the demographics and evolutionary paths of ob- 
served QSO/AGN. 

The search for the local QSO relics via the study of their dy- 
namical influence on the surrounding stars and gas carried out since 
the launch of the Hubble Space Telescope (see Richstone et al. 
1998; Ferrarese et al. 2008, and reference therein) led ultimately to 
a shift of paradigm in our understanding of the role of black holes in 
the formation of structures in the Universe. The discovery of tight 
scaling relations between SMBH masses and properties of the host 
galaxies' bulges l|Magorrian et al. 1998] IGebhardt et al. 20001 
IFerrarese & Merritt 2000l ITremaine et al. 20021 
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IMarconi & Hunt 20031 [Hiiring & Rix"2004) [Hopkins et al. 2007a^ 
clearly pointed to an early co-eval stage of SMBH and galaxy 
growth. 

However, local scaling relations have proved themselves 
unable to uniquely determine the physical nature of the SMBH- 
galaxy coupling, and a large number of feedback models have been 
proposed which can reasonably well reproduce these relations 
(Silk &Rees 1998, Fabian 1999: Kauffmann & Haehnelt 2000, 
[Umemura 200i_; C avaliere & Vittorini 2002 : Wyithe & Loeb 2003, 

[Granato e t al. 20041 |Di Matteo, Springel & Hernquist 2005 

[King 2005| | Sazonov et al. 20051 IEscala 2006t . Most likelv. it will 
only be by studying in detail the redshift evolution of these scaling 
relations that we will gain critical insight on the nature of the 
interaction between growing black holes and their host galaxies 
and on the nature of the accretion process and of jet formation 
in various regimes. In fact, many have recently attempted to 
push the study of scaling relations to 2 > 0, despite the large 
observational biases ^erloni, Rudnick & Di Matteo 2004; 
IMcLure et al. 20061 IShieTds et al. 2006. ,Peng et al. 2005| 
[Wu2007; ITreu et al. 20071 ISalviander et al. 20071 

rShen et al. 20081 |Schranun, Wisotzki & Jahnke 2008) , but the 
question of if, and, in case, by how much did these relations evolve 
with redshift has so far eluded any unambiguous answer. 

For all these reasons, it seems to us necessary to make the best 
effort possible to reconstruct the history of SMBH accretion in or- 
der to follow closely the evolution of the black hole mass function. 
That will serve as a benchmark to test various models for SMBH 
cosmological growth as well as those for the black hole-galaxy co- 
evolution. 

As opposed to the case of galaxies, where the direct relation- 
ship between the evolving mass functions of the various morpho- 
logical types and the distribution of star forming galaxies is not 
straightforward due to the never-ending morphological and pho- 
tometric transformation of the different populations, the situation 
is much simpler in the case of SMBH. For the latter we can as- 
sume their evolution is governed by a simple continuity equa- 
tion jCavaliere, Morrison & Wood 1971||S"mall & Blandford 19921 
lYu & Tremaine 2002t . where the mass function of SMBH at any 
given time can be used to predict that at any other time, provided 
the distribution of accretion rates as a function of black hole mass 
is known (see section [5] for more details). An analogous continu- 
ity equation could in principle be solved for the evolution of black 
hole spin (and thus for the overall accretion efficiency), whose time 
derivative, however, depends both on the accretion rate and on the 
merger rate evolution (Berti and Volonteri 2008; King, Pringle and 
Hofmann 2008). 

Consequently, the simple physical and mathematical structure 
of a black hole implies that their intrinsic properties are essentially 
determined by three physical quantities only: mass, spin and accre- 
tion rate. Leaving aside a few open issues related to the physics 
of strongly super-Eddington flows, accretion theory is a mature 
enough field to allow a relatively robust mapping of observables 
(spectral energy distributions, variability patterns, etc.) into these 
physical quantities associated to SMBH growth. This reduces the 
problem of unveiling the census of growing black holes to that of 
a wavelength-dependent assessment of selection biases, bolometric 
corrections, and accurate determinations of the distribution of ex- 
trinsic observational properties, the level of intervening obscuration 
being the main one. 

The fact that obscuration is a crucial factor in de- 
termining the observed properties of AGN was rec- 
ognized early after the discovery of the X-ray back- 



ground (XRB) radiation ( |Setti & Woltjer 1989) . Sub- 
sequent generation of synthesis models of the XRB 
(Madau, GMsellini & Fabian 1993[ [Comastri et a l. 19951 

Oilman & Fabian 19991 [Gilli, Salvati^fc Hasinger 2001[ 

[Gilli, Comastri & Hasinger 2007[ l, following closely on the 
footsteps of increasingly larger and deeper surveys, have progres- 
sively reduced the uncertainties on the absorbing column density 
distribution, which, coupled to the observed X-ray luminosity func- 
tions, provide now an almost complete census of the Compton thin 
AGN (i.e. those obscured by columns A^h < o-^^ ~ 1.5 x 10^'* 
cm^^, where ctt is the Thomson cross section). Indeed, this 
class of objects dominates the counts in the X-ray energy band 
where almost the entire XRB radiation has been resolved into 
individual sources jHasinger et al. I998[ IRosati et al. 20021 
[Alexander et al. 20031 IMoretti et al. 20031 [^rsley et al. 2005l l. 
Accurate determinations of the XRB intensity and spectral shape 
(for the most recent results see Churazov et al. 2007; Ajello et al. 
2008), coupled with the resolution of this radiation into individual 
sources ( [Brandt & Hasinger 2005[ >, allow very sensitive tests of 
how the AGN luminosity and obscuration evolve with redshift (the 
so-called 'X-ray background synthesis models' see e.g. Gilli et 
al. 2007, and references therein). In much the same way, accurate 
determinations of the local SMBH mass density and of the AGN 
(bolometric) luminosity functions, coupled with accretion models 
that specify how the observed AGN radiation translates into a 
black hole growth rate, allow sensitive tests of how the SMBH 
population (its mass function) evolves with redshift. This is the 
aim of this work, which, by analogy, we name 'AGN synthesis 
modelling' . 

Reconstructing in detail the history of accretion onto super- 
massive black holes is also of paramount importance for under- 
standing the effects that the energy released by growing black holes 
has on its environment. As we mentioned before, the various re- 
lations between black hole mass and the physical properties of 
host galaxy's bulge suggest that the growth of black holes and 
large scale structure is intimately linked. Moreover, X-ray obser- 
vations of galaxy clusters directly show that black holes deposit 
large amounts of energy into their environment in response to radia- 
tive losses of the cluster gas JBirzan et al. 2004TlFabian et al. 20061 
Allen et al. 20061. From these (and other) studies of the cavities, 
bubbles and weak shocks generated by the radio emitting jets in the 
intra-cluster medium (ICM) it appears that AGN are energetically 
able to balance radiative losses from the ICM in the majority of 
cases ( [Rafferty et al. 2006) . Also, numerical simulations of AGN- 
induced feedback have recently shown that mechanical feedback 
from black holes may be responsible for halting star formation in 
massive ellipticals, perhaps explaining the bimodality in the color 
distribution of local galaxies ( Spring eTet al. 2005|l , as well as the 
size of the most massive ellipticals ( De Lucia & B laizot 2007 ). 

All these arguments, however, hinge on the unknown effi- 
ciency with which growing black holes convert accreted rest mass 
into kinetic (ekin) and/or radiative (trad) power. Constraints on 
these efficiency factors are therefore vital for all models of black 
hole feedback, and we aim here at pinning them down as accurately 
as currently possible. 

The structure of this paper is the following: we will first start 
with a description of our current knowledge about the various 
modes in which a black hole may grow (§[2j. This will be used to 
introduce the function that relates the observed (bolometric) lumi- 
nosity of an AGN to the true accretion rate. In section[3]we will out- 
line the methods used to calculate the evolution of the SMBH mass 
(and accretion rate) function, while in section [4l we will spell out 
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and discuss the fundamental assumptions made in this work, in par- 
ticular those about the local mass function of SMBH, the adopted 
bolometric correction and the incidence of obscured AGN. Our re- 
sults on the growth of the SMBH population will be presented in 
§|5] There, we will first discuss "integral" constraints 15.1b and 
then "differential" ones (§ |5.2|5.4| l, which will reveal interesting 
features of the so-called anti-hierarchical (or 'downsizing') evolu- 
tion of SMBH. Section |6] will instead be devoted to the study of 
the kinetic luminosity function of AGN and of its evolution. At the 
end of it (§|6j2} we will be able to assess the relative importance of 
kinetic and radiative feedback from growing black holes as a func- 
tion of both mass and redshift. Finally, we will discuss our results 
and draw our conclusions in section[7] 

Throughout the paper we will use concordance cosmological 
parameters of Q,m = 0.3, Q,k ~ 0.7, Ho — 70kms^^ Mpc^^. 



2 THE MODES OF BLACK HOLE GROWTH 

Galactic (stellar mass) black holes in X-ray binaries, either of 
transient or persistent nature, are commonly observed undergo- 
ing so-called transitions, i.e. dramatic changes in their spec- 
tral and variability properties associated to changes in the way 
matter accretes onto them. The study of these various spectral 
states, either in individual objects or in ensemble of sources, 
has solidified in the last few years into a general framework 
of the various physical modes in which accretion onto a black 
hole can take place (Remillard & M cClintock 2006l|Malzac 20071 
[Done, Gierliiiski fe'Kubota 2007) . We discuss these modes, and 
their AGN analogues, in this section, paying particular attention 
to their partition in terms of the value of the Eddington ratio 
A = I/boi/I/Edd, where Lboi is the bolometric luminosity and 
LEdd = 47rGA/BHmpc/crT ~ 1.3 x 1O^®(Mbh/M0) erg s^^ 
is the Eddington one. 

There are at least three well defined spectral states. In 
the low/hard state the emission is dominated by a hard X-ray 
power-law with an exponential cutoff at about few 100 keV. Ra- 
dio emission is always detected in this state ( IFender 20061 1 usu- 
ally with flat or inverted spectrum, which is commonly inter- 
preted as due to a compact, persistent, self-absorbed jet of low 
intrinsic power. Interestingly, the radio luminosity of these jet 
cores seems to correlate quite tightly with the X-ray luminos- 
ity jGallo, Fender, & Pooley 2003^ . The spectrum of the higlVsoft 
state, instead, is dominated by a thermal component likely origi- 
nated in a standard Shakura & Sunyaev (1973) accretion disc, while 
in the intermediate (or very high) state, associated with transitions 
between hard and soft states, and often occurring at a source's high- 
est flux level, both a thermal and a steep power-law component 
substantially contribute to the spectrunjj. Maccarone (2003) has 
shown that the state change in these systems generally occurs at 
(bolometric) luminosities of about A ~ a few x 10^^. The hard-to- 
soft transition is also accompanied by a "quenching" of the steady 
radio emission observed in the low/hard state, following rapid flar- 
ing events in the radio band, often associated with fast ejection of 
matter moving at relativistic speeds and emitting optically thin syn- 
chrotron radiation ( [Fender, Belloni & Gallo 2004^ . 

Given the many analogies between the high-energy spectra 

^ In fact, the overall evolution of transient black holes does show clear ev- 
idences of hysteresis ( [Miyamoto et al. 1995[ IBelloni et al. 2005t . with the 
intermediate spectral states having different luminosity when sources enter 
or leave the outburst. We neglect this further complication here. 



of galactic black holes and AGN, a similar phenomenology has 
long been searched for in active supermassive black holes. Mer- 
loni, Heinz and Di Matteo (2003; MHD03), following earlier sug- 
gestions by Falcke and Biermann (1996), presented a thorough 
discussion on the scale-invariant properties of black hole coupled 
accretion/jet systems, based essentially on a generalization of the 
radio-X-ray correlation of low/hard state X-ray binaries to SMBH 
of varying masses. This led to the discovery of a so-called "funda- 
mental plane" of active black holes (see also Falcke, Kording and 
Markoff 2004; Kording, Falcke and Corbel 2006; Merloni et al. 
2006; Fender et al. 2007). 

This correspondence between low/hard state X-ray binaries 
and low-luminosity AGN is confirmed by a number of recent 
studies. Detailed multi-wavelength observations of nearby galaxies 
have revealed a clear tendency for lower luminosity AGN to be 
more radio loud as the Eddington-scaled accretion rate decreases 
(see Ho 2002, and references therein; Nagar, Falcke and Wilson 
2005). The above mentioned scaling relations with black hole 
X-ray binaries have thus helped to identify AGN analogues 
of low/hard state sources, in which the radiative efficiency 
of the accretion flow must be low (Churazov et al. 2005] 

Pellegrini 2005a[ |Chiaberge, C^etti & Macchetto 2005] 

Hardcastle, Evans & Croston 2007^ . 

In this work, we will use the kinetic jet power as a key phys- 
ical variable that identifies and characterizes the nature of various 
accretion modes. In fact, Merloni and Heinz (2007), using a sample 
of nearby AGN for which BH masses, nuclear radiative and kinetic 
luminosities and outer boundary condition (i.e. accretion rate at the 
Bondi radius) were simultaneously known, have found a clear re- 
lationship between Eddington-scaled kinetic power and bolometric 
luminosity, given by: log(Vykin/i'Edd) = (0.49 ± 0.07) log A - 
(0.78 ± 0.36). The measured slope suggests that these objects are 
in a radiatively inefficient accretion mode. Moreover, the observed 
correlations are in very good agreement with theoretical predic- 
tions of adiabatic accretion models with strong outflows, the rel- 
ative strength of which increases with decreasing accretion rate 
( [Blandford & Begelman 1999i|BIandford & Begelman 2004[ >. 

This is a very important finding, as the existence of jet 
dominated accreting SMBH bears important consequences for 
our understanding of the AGN feedback in a cosmological con- 
text. The specific importance of mechanical (as opposed to ra- 
diative) feedback for the heating of baryons within the deep- 
est dark matter potential wells has recently been acknowledged 
by semi-analytic modelers of cosmological structure formation 
dCroton et al. 20061 IBower et al. 20 06). Within these schemes, be- 
cause the bright quasar population peaks at too early times, a so- 
called "radio mode" of SMBH growth is invoked in order to regu- 
late both cooling flows in galaxy clusters and the observed sizes 
and colors of the most massive galaxies in the local Universe 
( [Springel et al. 2005[ICroton et al. 200"6t . 

Here, we also consider such a mode of accretion, where most 
of the released energy is in kinetic form, but we uniquely as- 
sociate it to black holes accreting below a certain critical rate 
in Eddington units, rather than to the state of the fuelling gas 
( [Hardcastle, Evans & Croston 2007^ , or to the properties of the 
very large scale environment (as, for example, the host dark mat- 
ter halo size). In order to avoid confusion, we call this "low kinetic 
(LK)" mode, to distinguish it not only from the bright "high ra- 
diative (HR)" mode (the 'quasar' mode of the recent cosmology 
jargon, associated to radio quiet QSOs, analogous to the high/soft 
state of X-ray binaries), but also from the bright radio loud quasars 
(of which 3C273 is the prototype), characterized by both radiatively 
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efficient accretion flows and powerful jets, that we term "high ki- 
netic (HK)" mode, as we discuss below. 

The most luminous sources, as those falling into the stan- 
dard definition of Quasars and broad lined AGN, are most 
likely accreting at a high rate, close to the Eddington limit 
( |McLure & Dunlop 2004{ IKollmeier et al. 20061 1. Also their spec- 
tral energy distribution, usually dominated by the so-called Big 
Blue Bump (BBB: quasi thermal UV emission from an optically 
thick standard accretion disc, see e.g. Malkan 1983; Laor 1990), 
indicates that above a certain critical value of the ratio Acr, an accre- 
tion mode transition should take place to what is usually described 
as a standard, geometrically thin and optically thick accretion disc 
(Shakura & Sunyaev 1973). Evidence for such a transition taking 
place in AGN has been presented by Maccarone, Gallo & Fender 
(2003), Jester (2005) and Greene, Ho and Ulvestad (2006). 

Thus, according to our scheme, black holes accreting above 
the critical rate come into two different physical states, one with 
(HK), and one without powerful jets (HR). The reason for this di- 
chotomy (and even its reality, see e.g. Cirasuolo et al. 2003) has 
been subject to countless arguments since the very early years of 
the first radio loud QSOs discoveries, but it is still matter of debate 
(for example, see Sikora, Stawarz & Lasota 2007; Kaiser & Best 
2008; Blundell 2008 for some recent ideas), with the two most pop- 
ular options for the observed dichotomy being the influence of the 
black hole spin and a change in accretion mode. 

As a working hypothesis, in line with what originally proposed 
in MHD03 (see also Nipoti et al. 2005; Kording, Jester and Fender 
2006; Blundell 2008), we consider here that powerful radio jets are 
episodic and transients events during the time a black hole spends 
accreting above the critical rate (i.e. close to the Eddington limit). 
This idea, that is easy to incorporate into any population synthesis 
scheme, bears resemblance with what actually observed for tran- 
sients black hole X-ray binaries, which often displays very power- 
ful and rapid radio flares at the peak of their outbursts (see Fender 
et al. 2004, and references therein). 

2.1 A schematic view of radiative, liinetic and accretion 
efficiency 

We adopt a simple, physically motivated functional form for the A 
vs. m function (where m is the accretion rate onto the black hole 
in Eddington units m = rjMc^ / L-^dd, vvith rj accretion efficiency) 
of a broken power-law, bridging the low accretion rate (radiatively 
inefficient) regime (A cx m^) and the high accretion rate one. For 
radiatively efficient sources, we should assume, by definition, that 
the bolometric luminosity is simply proportional to the accretion 
rate, A oc m. The overall normalization is found by imposing conti- 
nuity at Acr and depends on the adopted value of this critical accre- 
tion rate. Based on the results of Maccarone et al. (2003), Greene 
et al. (2006), Merloni and Heinz (2007), we adopt here a value of 
Acr = 3 X 10~^ (see also the discussion in Hopkins, Narayan & 
Hernquist 2006). A schematic view of these various modes in the 
accretion rate vs. released power plane is shown in Fig.[T] The Fig- 
ure represents qualitatively the main features of the three-modes 
scheme we have outlined here: below the critical rate, only one 
mode (LK) is allowed, where most of the power emerges in kinetic 
form and there is a quadratic relationship between accretion rate 
and (bolometric) radiative luminosity; above it, two different modes 
are possible: one where kinetic and radiative power are compara- 
ble and high (HK), and one where kinetic power is quenched (HR). 
The exact slopes of the kinetic power vs. accretion rate relation are 
set, implicitly, by the relationship between radio core luminosity. 
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Figure 1. A schematic view of the relationship between accretion rate and 
released power, both in units of the Eddington luminosity (an "accretion 
mode map"). Blue solid lines show the radiated power, red dashed ones 
the released kinetic power. The black horizontal dotted line mark the criti- 
cal Eddington ratio where a change of accretion mode is assumed to occur 
Below this rate, only one mode (LK) is allowed, where most of the power 
emerges in kinetic form; above it, two different modes are possible: one 
where kinetic and radiative power are comparable and high (HK), and one 
where kinetic power is quenched (HR). The lower panel shows the corre- 
sponding ratios between kinetic and radiative power as a function of accre- 
tion rate. 



kinetic power and X-ray luminosity (the "fundamental plane" of 
active black holes), as well as by the hard X-ray bolometric correc- 
tion function, that we discuss in §[3]and|4] 

This simple relation, together with a bolometric correction, al- 
lows the determination of the dimensionless accretion rate for any 
object for which both the nuclear luminosity, Li, (in any band) and 
the black hole mass are measured. Conversely, the accretion lumi- 
nosity of a black hole of mass M, growing at a rate M can be 
obtained from the above relationship. In § [3] we will describe in 
detail the method we employ to unveil the evolution of the black 
hole mass function based on the evolution of the intrinsic hard 
X-ray (2-10 keV) luminosity function, making use of the relation 
Li = Li{Ai, Al) discussed above. 

As a final remark in this section, we come back here to the 
well known distinction between accretion efficiency, and radiative 
efficiency, that will be important in the following. The former repre- 
sents the maximal amount of potential energy that can be extracted, 
per unit rest mass energy, from matter accreting onto the black 
hole. This quantity, 77(a), depends on the inner boundary condi- 
tion of the accretion flow only, and, in the classical no-torque case 
dNovikov & Thome 19731 1 is a function of BH spin, a, only, ranging 
from 77(0 = 0) ~ 0.057 for Schwarzschild (non-spinning) black 
holes to ri{a = 1) ~ 0.42 for maximally rotating Kerr black holes. 

On the other hand, the radiative efficiency, trad = L-hoi/Mc^ 
depends both on the accretion efficiency (i.e. on the inner boundary 
condition of the accretion flow) and on the nature of the accretion 
flow itself. In the case of a simple broken power-law shape that we 
have assumed for the A — m relation, we have that: 
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^rad 



; ?7/(rh) = 77 x 



1, m > rhc 

(rh/rhcr), rh < rhc 



(1) 



where rhcr = Acr is the critical Eddington- scaled accretion rate 
above which the disc becomes radiatively efficient. As we will 
show in § [5] the classical argument due to Soltan (1982) can 
be used to put constraints on the (mass weighted) average value 
of (erad), by comparing the evolution of AGN/QSO luminosity 
function with estimates of the local relic black hole mass func- 
tions. Constraints on {r](a)), and thus on the average spin of the 
SMBH population, can then be derived from the use of eq. (TJ 
( [Hopkins, Narayan & Hernquist 2006] l. The relationship between 
these average values (77) and (e^ad) depends both on the value of 
Acr (higher values of the critical accretion rate increase the number 
of objects which are in a radiatively inefficient regime, thus reduc- 
ing the ratio between radiative and accretion efficiencies, and vicev- 
ersa) and on the shape and evolution of the luminosity functions 
themselves. However, in the remain of the paper we will assume 
that eq. ([TJ holds, and will not discuss explicitly the dependence 
of the (r7)/{erad) ratio on Acr. Moreover, for the sake of clarity, 
we will refer in the following to the mass weighted averages of the 
efficiencies simply with the symbols 77 and e^ad- 



3 METHODS 

In the simplest case of purely accretion driven evolution (i.e. as- 
suming mergers do not play an important role in shaping the black 
holes mass function, but see § [7] for more on this point), the si- 
multaneous knowledge of the black hole mass function and of the 
Eddington ratio distribution function, '^m^\{z) (i.e. the number of 
objects with given mass M and Eddington ratio A per comoving 
volume and logarithmic interval), would completely determine the 
evolutionary solution for the mass function, as the two distribu- 
tions must be coupled via a continuity equation (Small & Blandford 
1992; Yu & Tremaine 2002; Marconi et al. 2004; Merloni 2004, 
hereafter M04; Shankar, Weinberg & Miralda-Escude 2008) 



9$M(M,t) 9[<&M(M,t)-(M(M,f))] 



dt 



dM 



(2) 



to be solved given an appropriate boundary condition. Here <1>m is 
the black hole mass function, and (M(M, t)) represents the aver- 
age accretion rate (in physical units) for a SMBH of mass M at the 
time t, and can be uniquely determined once the black hole Edding- 
ton ratio distribution function, <1?m,a(z), is known. 

The method used to solve eq. (|2j follows closely that applied 
in M04, with some modifications and improvements, as we de- 
scribe it in detail here. As in M04, we argue that the most robust 
constraint on the SMBH mass function at any redshift can be ob- 
tained by integrating the continuity equation backwards in time, 
using the local mass function as a boundary condition. Moreover, 
as opposed to all other previous calculation of the SMBH mass 
function evolution, we do not make any simplifying assumption 
about the distribution of the accretion rates for different black hole 
masses and redshifts, but we do calculate these distributions from 
the multi-wavelength luminosity functions of the AGN population. 
The guiding principle is that, within each of the three physically 
different modes of accretion described in the previous section (LK, 
HR, HK), there exist well defined relationships between an object 
luminosity in various bands and the physical parameters needed to 
solve eq. i.e. the black hole mass and its accretion rate. In sec- 
tion[2]we have already presented the relationship between accretion 
rate, mass and bolometric luminosity in these various modes. The 



second relationship, needed in order to determine the evolution of 
the AGN jet (and thus of their kinetic energy feedback) is that be- 
tween radio luminosity, bolometric luminosity and BH mass, that 
we discuss here. 

In the 'low kinetic' (LK) mode, i.e. for all sources with Ed- 
dington ratios A less than the critical value Acr, we adopt the 
now well established result, first described in MHD03, that a 
'fundamental plane' relation between (core radio) jet luminosity, 
(accretion-powered) X-ray luminosity and black hole mass can in- 
deed be defined for those objects (but see, for example, Panessa et 
al. 2007 for a contrasting view). We write it here as 



log Lb. = Crx log Lx + Crm log M + /3lk, 



(3) 



where Al is the black hole mass in solar units, Lr is the intrin- 
sic (i.e. un-beamed) radio luminosity of the jet core at 5 GHz and 
Lx stands for the intrinsic (i.e. unabsorbed) X-ray luminosity in 
the 2-10 keV band (both in units of erg s^^), which we uniquely 
relate to the bolometric luminosity by the relation of eq. (21) of 
Marconi et al. (2004). We adopt for the correlation coefficients the 
values recently determined in Merloni and Heinz (2007) based on 
a more accurate study of the effects of relativistic beaming, as well 
as a more accurate selection of the parent sample, as first discussed 
in Kording, Falcke and Corbel (2006). Specifically, we assume 
Crx = 0.62, 5rm ~ 0.55 and /3lk = 8.6. We also assume that this 
relationship has an intrinsic scatter of 0.6 dex, a value smaller than 
that found in MHD03, but in line with that discussed in Kording et 
al. (2006) for samples of jet dominated objects only. 

Let us now discuss the relationship between radio, X-ray lu- 
minosity and black hole mass for objects accreting at high rates. Ir- 
respective of the physical explanation preferred for the existence of 
a radio loud/quiet (HK/HR) dichotomy, the problem we face for our 
method is that there are no well defined relationships to be found in 
the literature between radio luminosity, mass and bolometric lumi- 
nosity which are specific for black holes accreting above the critical 
rate, Acr. One reason for this, as we have argued before in Merloni 
et al. (2005), is that, because such modes can exist only in very 
limited range of the parameter space (i.e. for Acr ^ A < 1, where 
A = 1 corresponds to the Eddington limit), it is almost impossible 
to unravel the mass from the accretion rate dependence, given the 
usual uncertainties on SMBH mass estimates and the likely pres- 
ence of intrinsic scatter in any such relation. Thus, our treatment of 
the HK sources (classical radio loud QSOs) should necessarily be 
considered at best as an educated guess. For the sake of simplicity, 
we assume that their intrinsic radio core luminosity scales with the 
bolometric luminosity, so we fix 



log Lr = log Lbol + Chm log M + /3hk,hr. 



(4) 



We work out the mass scaling ^hm (the same for HK and HR 
modes) and the normalization of the HK (radio loud) mode assum- 
ing continuity with the fundamental plane relation at the critical 
accretion rate, while for the HR (radio quiet) mode, we simply 
take a luminosity 10'^ times smaller (/3hr = /3hk — 3), so that 
the radio luminosity in the HR (radio quiet) mode is effectively 
quenched. We obtain ^hm = 0.02 (mass dependency is negligible), 
and /9hk = —4.9, in relatively good agreement with the standard 
radio loud SED of Elvis et al. (1994). Also for this relationship we 
assume an intrinsic scatter of 0.6 dex. 

Once fixed the relationship between Lr, Lx and M for all val- 
ues of the accretion rate (i.e. for the three main modes LK, HR and 
HK), we can proceed in extracting the distribution of accretion rate 
as a function of SMBH mass and redshift. In order to do so, we take 
the observed luminosity functions for Lr (intrinsic radio core emis- 
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sion) and Lx (intrinsic 2-10 keV luminosity) and the SMBH mass 
function. We then look for the joint distribution 'I'c(^h, Lx, M), 
i.e. the function that quantifies the differential number of objects 
(per unit comoving volume) with given mass AI, radio luminosity 
Lr and X-ray luminosity Lx- We do so by imposing the relations 
1(3) and l|4j, with their own intrinsic scatter, and minimizing the 
differences between the projections of ^'c(-f'R, Lx, M) onto the 
three axes identified by Lr, Lx and M and the observed luminos- 
ity and mass functions. For high accretion rate objects (i.e. those 
with A > Acr), we assume that the HK (radio loud) objects repre- 
sent only 10% of the total, so that the number ratio of HR/HK AGN 
is about 9/1. We also impose the Eddington limit by assuming that 
the number of objects with A > 1 are exponentially suppressed. By 
construction, then, our joint distribution function vE'c(Lr,, Lx, M) 
reproduces the observed luminosity (in the radio and X-ray bands) 
and mass functions, and is consistent with the scaling relations l[3) 
and Q. 

Once the joint distribution "ifc is found, it is straightforward 
to infer the accretion rate distribution function: each combination 
(Lr, Lx and M) correspond to a unique value of the accretion 
rate Af , and the number of objects per unit comoving volume (per 
logarithmic interval) with a given value of M and M is given by: 

*M j\/ = ^ 

dlogMdlogM 

= I *c[LR,Lx(M,M),M]^iHii2£dlogLR (5) 
J d\ogM 

where Lx = Lx{M, M) is the relation between intrinsic X-ray lu- 
minosity and mass and accretion rate discussed in section[2] Thus, 
the knowledge of the joint distribution function '^c, complemented 
with the theoretical relationship between luminosity and accretion 
rate for the various modes (X{M , M), see § |2]l is sufficient to 
derive the distribution of accretion rates for black holes of differ- 
ent masses. The average {M{M)) is what is then needed to solve 
eq. l|2j. 

In this way, the combined knowledge of the SMBH mass func- 
tion and of the AGN luminosity function in at least one band at a 
given moment in time completely determine the boundary condi- 
tions needed to integrate the continuity equation. In practice, start- 
ing from the mass function and the accretion rate distribution at a 
given redshift Zi, it is possible to derive the new black hole mass 
function at redshift Zi -\- dz, "1>m(M, Zi -I- dz), by just subtract- 
ing the mass accreted in the time interval dt — dz{dt/dz) calcu- 
lated according to the accretion rate function of redshift Zi. This 
new mass function can then be used together with the radio and 
X-ray luminosity functions at the same redshift, 0r(Lr, z + dz) 
and (j)x{Lx, z + dz), to obtain the new joint distribution function 
^'c(Lr, Lx, M, z+dz), and therefore the new accretion rate func- 
tion, and so on, until a solution of eq. l|2j is found in the interval 
Zi + Az, where the value of Az depends solely on the availability 
of good enough data (luminosity functions) to meaningfully deter- 
mine the joint distribution 'i'c- 



4 ASSUMPTIONS AND OBSERVABLES : M ASS AND 
LUMINOSITY FUNCTIONS 

According to the procedure we have outlined in §[3]above, in order 
to obtain meaningful constraints on the accretion rate distribution, 
we have to make the following assumptions: 



• The local BH mass function is known accurately down to a 
certain SMBH mass; 

• The bolometric corrections are known at all redshifts for AGN 
of all luminosities, so that meaningful bolometric luminosity func- 
tions can be derived; 

• The selection function of hard X-ray surveys provides us with 
a complete census of the AGN population, i.e. is not strongly bi- 
ased against obscured sources. Although this is obviously not the 
case for the so-called Compton Thick AGN, i.e. those obscured by 
such a large column (A*'h > cr^^) that also hard X-rays cannot be 
transmitted to the observer, we will proceed under the assumption 
that the Compton Thick AGN fraction can be meaningfully con- 
strained on the basis of the X-ray background models, and verify, a 
posteriori that their contribution to the overall growth is small; 

• The bolometric LF of AGN, once extrapolated to low enough 
luminosities, recovers the entirety of the local SMBH population, 
as described by the adopted mass function. This is equivalent to 
the statement that the AGN fraction approaches unity towards very 
low nuclear luminosities, once incompleteness of any specific LF 
is accounted for. 

We will discuss these issues in turn below. 



4.1 Local black hole mass function 

As it is well known, uncertainty in the local mass function 
determination propagates linearly in the computation of the 
average accretion efficiency in Soltan-type of arguments. This 
issue has been dealt with extensively in the literature fSoltan 1982) 
IMarconi et al. 20041 [Merloni, Rudnick & Di Matteo 2004). 

Recently, progress have been made towards a more secure com- 
putation of the local mass function of supermassive black holes, 
based on the various scaling relations. For example, Graham and 
Driver (2007) have shown how most previous inconsistencies 
among local black hole mass densities, pbh,o, computed using 
different scaling relations were due to incorrect accountancy 
for the Hubble constant dependence. Here, we will not dwell 
into a lengthy discussion about the available constraints on the 
local distribution of SMBH masses, for which we refer the 
reader to the comprehensive treatment presented in Shankar et 
al. (2008b), who have carried out a detailed comparison among 
SMBH mass functions derived with different methods. Rather, 
we briefly describe in the following our choice of the initial 
condition for eq. Although uncertainties do remain, with the 
local black hole mass density value estimated within the range 
pBH.o = (3.2 - 5.4) X 10^ M0 Mpc"^ (for JLq = 70 km 
s~^ Mpc^'^), the agreement between various measurements is 
encouraging. We will adopt a specific analytic expression for the 
SMBH mass function, coincident with the central value within 
the uncertainty range in Shankar et al. (2008b). Our adopted mass 
function is computed as the convolution of a Schechter function 
5'm with a Gaussian scatter of 0.3 dex, to account for the intrinsic 
scatter of the scaling relations from which the mass function 
is derived. For the Schechter function, we adopt the following 
parametrization 

with parameters: 0, = 10"^, log M* = 8.4 and a = -1.19. With 
this choice, the local black hole mass density is pbh,o ~ 4.3 x 10'' 
M0 Mpc~^ and the total number density of black holes above 5 x 
10^ solar masses is about nen — 1-3 x 10^^ Mpc^^. 
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Figure 2. Redshift evolution of the accretion rate density onto supermassive 
black holes 'I'bharX^), multiplied by a constant (equal to 10'') is plotted 
as a shaded area alongside the best fit to a large compilation of data for the 
star formation rate (SFR) density (blue dashed, from Hopkins and Beacom 
2006). The shaded area represents the uncertainty in the black hole accre- 
tion rate density deriving from the obseiTational uncertainty on the AGN 
luminosity functions. 



4.2 Bolometric corrections and luminosity functions 

The cosmic evolution of the SMBH accretion rate and its associ- 
ated mass density can be calculated from the bolometric luminosity 
function of AGN: <^(Lboi,z), where Lboi = ^ra.dM(? is the (in- 
trinsic, i.e. before absorption) bolometric luminosity produced by 
a SMBH accreting at a rate of M with a radiative efficiency trad • 
In practice, the accreting black hole population is always selected 
through observations in specific wavebands. Crucial is therefore the 
knowledge of two factors: the completeness of any specific AGN 
survey, and the bolometric correction needed in order to estimate 
Lboi from the observed luminosity in any specific band. 

As for the accuracy of our knowledge of the bolometric 
correction, we refer the reader to the studies of Marconi et 
al. (2004); Richards et al. (2006); Hopkins et al. (2007b) and 
Shankar et al. (2008b). All of them consistently demonstrate 
that a luminosity dependent bolometric correction is required in 
order to match type I (unabsorbed) AGN luminosity functions 
obtained by selecting objects in different bandfl No signifi- 
cant trend with redshift is instead evident, similarly to what is 
known for the spectral properties of type I AGN in the X-ray 
band alone ([Vignali, Brandt & Schneider 2003] ITozzi et al. 20061 
ISteffen et al. 2006t . Here we will adopt, as fiducial inputs, the bolo- 
metric corrections of Marconi et al. (2004) and the latest observed 
2-10 keV luminosity function of Silverman et al. (2008), which pro- 
vides the best available constraint on X-ray selected AGN at high 
redshift. In order to account for the measurement uncertainties in 



^ This conclusion, however, has been recently challenged by Vasudevan & 
Fabian (2007), who ague that the 2-10 keV X-ray bolometric correction is 
better correlated with Eddington ratio rather than with luminosity. 



a simple and straightforward way, we will use two different ana- 
lytical parametrization of the observed 2-10 keV luminosity func- 
tion, as discussed in Silverman et al. (2008): a luminosity depen- 
dent density evolution model (LDDE; Ueda et al. 2003; Hasinger et 
al. 2005) and a modified-pure luminosity evolution (MPLE; Hop- 
kins et al. 2007b). The analytic expressions for these models can be 
found in eqs. (5-11) and in Table 4 of Silverman et al. (2008). 

4.3 Accounting for obscured sources 

Hard X-ray selected samples of AGN represent the most unbiased 
census of accretion activity in the Universe, due to the limited effect 
that (Compton thin) gas obscuration has on the emerging luminos- 
ity in this band. However, observed 2-10 keV luminosities of AGN 
are still affected by absorption, and we need to quantify this effect 
when looking for the intrinsic bolometric luminosity. 

Assuming an average X-ray spectral model of the Gilli et al. 
(2007) X-ray background synthesis model, with fixed power-law 
slope of r = 1.9, we can compute the ratio of emerging to in- 
trinsic 2-10 keV luminosity, g = I/x.obs/ix, as a function of 
the absorbing column density of cold gas, A^h (in cm~^). Using 
XSPEC (Arnaud 1996) we have computed the effect of various ob- 
scuring columns on the emerging spectrum and then fitted the ob- 
served relation between g and A'h with a fourth order polynomial in 
A/h = log A^H — 21, and found the following analytical expression, 
accurate to better than 20% in the whole range 10^^ < A^h < 10^'*: 

g{Afn) ~ 0.99+0.018A^h-0.059A/'h^-0.1A/'h^ + 0.028A/'h''(7) 

For Compton thick objects (log A''h > 24), we have simply 
assumed that only a scattered component of 0.02 times the intrinsic 
flux is detectable (see Gilli et al. 2007). 

Then, for any given observed 2-10 keV luminosity func- 
tion </>'(Lx), we can deconvolve the effect of obscuration if 
the distribution of column densities f{Nu,Lx) as a function 
of luininosity is known. It is now well established that the 
fraction of absorbed sources depends on the intrinsic X-ray 
luminosity of an AGN (Ueda et al. 20031 |La Franca et al. 20051 
[Treister & Urry 2005[[Barger et al. 2005] [Hasinger 2008| ). Follow- 
ing again the latest incarnation of the XRB synthesis model of Gilli 
et al. (2007), we can write the ratio of absorbed to unabsorbed 
sources, R{Lx) as 



1 + 2.7e" 



(8) 



with Lc = 10 , so that the ratio of obscured to unobscured 
sources tends to unity at high (QSO) luminosities and to 3.7 at low 
luminosities. 

The distribution of absorbing columns is taken to be an in- 
creasing function of the column density for Compton thin objects, 
and, for the lack of better observational constraints, flat for Comp- 
ton thick sources, once again fully consistent with the XRB synthe- 
sis model of Gilli et al. (2007). We thus define 



where 

r4.5 



fi(^x) 

o /I R{L^) 
■^^l-H(Lx) 

fixed by 



A4t<3 
AAh > 3 



(9) 



A is fixed by the overall normalization 

/;V(A^H,Lx)dVH = l. 

Finally, we obtain the following expression for the intrinsic 
2-10 keV AGN luminosity function as a function of the observed 
(un-corrected) one: 



<^(Lx) ~ (l)'{L^) 



84, 



91og Lx 



1 - <?(A/'H)],f(7VH,Lx)dArH(10) 
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and this is what we use in our calculation of 'i'c (section[3](. 

4.4 Do all SMBH have an "active" counterpart? 

Traditionally, the presence of an accreting supermassive black hole 
in the nucleus of a galaxy have been associated to a number of 
unambiguous signs of activity: high excitation broad and/or nar- 
row emission lines, hard X-ray luminosity exceeding 10^^ erg s~^, 
the presence of relativistic jets, etc. Within this view, it is com- 
monly held belief that the fraction of AGN in any randomly sam- 
pled population of galaxies is relatively small. This is in obvi- 
ous contrast with the idea that SMBH are ubiquitous in the nuclei 
of galaxies (at least in the nearby Universe), and have prompted 
many authors to search for the "trigger" that activate an other- 
wise quiescent black hole. However, it must be obvious that ev- 
ery observational definition of AGN, like those mentioned above, 
suffers from strong selection biases, in particular at low intrinsic 
luminosities. It is in fact well known that high-resolution, high- 
sensitivity studies of nearby galaxies not considered to be powerful 
AGN reveals that a very high fraction of them do indeed harbour 
SMBH showing clear signs of accretion. For example, Nagar et 
al. (2005) and Filho et al. (2006) have used high-resolution radio 
imaging of a magnitude limited sample of bright nearby galaxies 
(the Palomar sample) to show that at least a quarter (and possi- 
bly as much as half) of them are true AGN, albeit of very low 
luminosity. Similar results have been obtained by optical spec- 
troscopic studies jHo, Filippenko, & Sargent 1991\ and by sensi- 
tive Chandra X-ray observations dHo et al. 2001||Gallo et al. 20081 
I S antra, Sanders & Fabian 2007 [ i. 

Furthermore, we notice that at z — 0, the low-luminosity 
slopes of both radio ( [Nagar, Falcke & Wilson 2005 [ > and X-ray 
dSilverman et al. 2008b luminosity functions of AGN are steeper 
than the low mass end slope of the SMBH mass function. There- 
fore, if we limit ourselves to the study of SMBH with (relic) masses 
larger than a given value, only a minor extrapolation of the observed 
luminosity functions towards low luminosity is sufficient to reach 
total number densities comparable to that of the local SMBH. 

For example, for our adopted luminosity functions, the num- 
ber of SMBH with M > 5 x 10^ Mq is equal to the num- 
ber of objects with intrinsic (i.e. un-absorbed) 2-10 keV luminos- 
ity larger than a few xlO*" erg s^^ and with intrinsic (i.e. de- 
beamed, see § I6.lt 5 GHz radio core luminosity greater that a 
few xlO^* erg s~^. From hard X-ray surveys of the local Uni- 
verse (Sazonov et al. 2007J, as well as from XRB synthesis models 
( [Gilli, Comastri & Hasinger 2007| l, we know that at these very low 
luminosities the number of X-ray AGN is dominated (in a ratio 
of approximately 4:1) by heavily obscured objects, so we expect 
that sensitive Chandra surveys reaching down to these luminosi- 
ties should unveil nuclear accreting black holes in a large fraction 
of nearby galaxies, as indeed found by Gallo et al. (2008). On the 
other hand, radio surveys are insensitive to dust obscuration (al- 
though affected by relativistic beaming), and indeed, based on the 
works of Nagar et al. (2005) and Filho et al. (2006) it is reasonable 
to expect that essentially all nearby galaxies will show compact, 
high brightness temperature cores with Lsghz > 10^^ erg s~^. 

In our view, this justifies the assumption, implicit in all our 
subsequent calculations, that there is no distinction between SMBH 
and AGN, in the sense that every supermassive black hole is active, 
at some level, like the SMBH at the center of our galaxy clearly 
demonstrates. Thus, the SMBH population coincides with that of 
all actively accreting black holes, once the broad distribution of 
accretion rates is accounted for . This is a very important point, at 



variance with most current analytic, semi-analytic and numerical 
models of SMBH growth that often assume a fixed Eddington ra- 
tio, and the consequences thereof will become manifest more than 
once in the remain of the paper. In particular, we will discuss the 
observational implications for the determination of AGN fraction 
in different surveys in section [54l 



5 UNVEILING THE GROWTH OF SUPERMASSIVE 
BLACK HOLES 

In this section, we concentrate on the growth of the supermas- 
sive black hole population. We will proceed from the most general 
to the more detailed results, namely from the study of integrated 
quantities (SMBH mass and accretion rate density, radiative energy 
density, average radiative efficiency, etc.) and their redshift evolu- 
tion, to differential ones, focusing on the full mass and accretion 
rate functions, and their evolution. The first part (§|5TTJ will shed 
light on some very general properties of the cosmological growth of 
SMBH, while the second (i) |5.2|5.4"l l will reveal interesting features 
of their so-called anti-hierarchical behaviour. 

In this section, as well as in the next one, all results will be 
shown accounting for the intrinsic uncertainties of the adopted lu- 
minosity functions. We estimate that these uncertainties can be 
evaluated by comparing different analytic parametrization of the 
same data sets (the LDDE and MPLE parametrization for the hard 
X-ray luminosity function of Silverman et al. 2008, 14.21 and the 
two alternative parametrizations for the flat-spectrum radio lumi- 
nosity function of Dunlop and Peacock (1990) and De Zotti et al. 
(2005), see §[6l(. Then, we will solve the continuity equation for the 
SMBH evolution according to the method described in section [3] 
for the four possible combinations of X-ray and radio luminosity 
functions, and present the results as shaded areas encompassing the 
entire range spanned by these four solutions. 

5.1 Integral constraints, mass density and radiative efficiency 

The redshift evolution of the black holes accretion rate (BHAR) 
density can be easily calculated from the intrinsic X-ray luminosity 
function l llOb as follows: 

*BHArW= r ' ^^ad)I.bol(Lx) ^^^^^ ^)rflogLx, (11) 
Jo Ei-adC^ 

where ix is the intrinsic X-ray luminosity in the rest-frame 2- 
10 keV band, and the bolometric correction function Lboi(ix) is 
given by eq. (21) of Marconi et al. (2004). 

In Figure|2]we show the redshift evolution of ^'bhar(2), mul- 
tiplied by a constant (equal to 10'^) in order to visually compare its 
shape with that of the star formation rate (SFR) density, which we 
take from the best fit to a large compilation of data from Hopkins 
and Beacom (2006). The evolution of the BHAR density of Fig. ([2J 
has been computed for erad — 0.07, but changing this value would 
simply result in a change of the overall normalization. 

By construction, the computed SMBH accretion rate density 
includes both absorbed and unabsorbed sources as well as Compton 
thick AGN (§ |43J, based on the prescription of Gilli et al. (2007) 
that uses, as almost unique observational constraint, the shape and 
normalization of the XRB radiation, and is thus somewhat degen- 
erate in the luminosity-redshift plane (see, however, Worsley et al. 
2005). As a sanity check, we can compute the total contribution 
of Compton thick sources to the total mass accumulated in SMBH 
today, which turns out to be limited to less than about 20%. This 
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Figure 3. Top panels: the redshift evolution of the SMBH mass density. 
Bottom panels, the evolution of the average Eddington ratio of the SMBH 
population, defined in eq. (13). On the left we show the results of a calcu- 
lation performed fixing the accretion efficiency to r; = 0.1, which, for the 
particular model of accretion modes we assume correspond to an average 
radiative efficiency of e^ad — 0.086, while on the right we show the case 
of 7? = 0.083, corresponding to Erad — 0.07. See text for details. 



confiims that our results are not too sensitive to the uncertainties in 
the luminosity and redshift distribution of Compton thick AGN. 

The integrated SMBH mass density as a function of redshift 
can then be calculated starting from a local value for pbh.o (see 
sectionl4Tt: 



Pbh(z) 
Pbh,o 



1 



BHAR 



{z') dt 



Pbh.o 



dz' 



dz' . 



(12) 



In this way, for any given <?!>(I/x, z) and bolometric correction (thus 
^bhar(2)), the exact shape of pbh(z) then depends only on two 
numbers: the local black holes mass density pbh.o and the (av- 
erage) radiative efficiency e^ad. This fact has an interesting, and 
somewhat surprising, consequence: even before making any as- 
sumption about the distribution of the Eddington ratio in the AGN 
population, we are able to constrain the (mass weighted) average 
radiative efficiency, in the following way. A firm lower limit to 
trad can be established using the Soltan argument. Specifically, it is 
enough to look for the value of the radiative efficiency that makes 
the integral -^-^f^^^ i^dz' of eq. ^ larger than unity. For 
our choice of the local SMBH mass function and AGN bolomet- 
ric luminosity function this lower limit is approximately equal to 
0.065. 

However, when starting from the local mass density to com- 
pute its evolution, the shape of pbh {z) turns out to be very sensitive 
to the average radiative efficiency. The very fact that, in much the 
same way as the SFR density, also the black hole accretion rate den- 
sity, \['bhar(z) has a well pronounced maximum (around z ~ 2), 
implies that, if the adopted value of trad were too high, the black 
hole mass density would not decrease any longer towards higher 
redshift. This is clearly shown in the upper two panels of Figure[3] 
where we show the evolution of the total SMBH mass density com- 



puted for two different values of the accretion efficiency (and, cor- 
respondingly, of the radiative efficiency): ?7 = 0.1 (corresponding 
to e^ad ^ 0.086) and 77 = 0.083 (e^ad 0.07). 

The same issue can be considered from a different point of 
view. Let us now define a global, mass weighted average Eddington 
ratio for the entire SMBH population as: 



K{z) = 4.38 X 10 Erad 



PBh(2) 



(13) 



The redshift evolution of A(z) is shown in the two bottom pan- 
els of Fig.[^ There it is apparent that, while for the low radiative 
efficiency case (lower right panel) the average Eddington ratio is al- 
lowed to increase up to the highest redshift for which reliable X-ray 
luminosity data are available (z ~ 4), increasing the efficiencies 
by less than 20% modifies the evolution of the A substantially (see 
bottom left panel), being it almost flat, and possibly decreasing at 
2 > 3, an obvious consequence of the flattening of the correspond- 
ing mass density (top left panel). The physical explanation for this 
is straightforward: if we calculate the evolution of the SMBH mass 
density with a too high radiative efficiency, we fail to account for 
the mass locked up in SMBH in the local Universe, and we are 
forced to explain a "primordial" population of very weakly active 
(or even dormant) black holes, of substantial density in place by 
2 ~ 4 (in the example of Fig. [3] left panel, amounting to about 
10^ Mq Mpc^'^). This puts already some strain on existing mod- 
els for the black holes seed population (see e.g. Trenti & Stiavelli 
2006; Volonteri, Lodato and Natarayan 2007). Models in which 
the primordial black holes originate from stellar mass progenitors 
(remnants of the first, PopIII, stars; Abel, Bryan & Norman 2000; 
Bromm, Coppi & Larson 2002), predict negligible BH mass den- 
sity at high redshift (say, less than lO"* Mq Mpc~^ at 2 > 5), 
and are incompatible with an average, mass weighted radiative ef- 
ficiency larger than about 0.075 (in very good agreement with the 
recent estimates of Shankar et al. 2008b), or, correspondingly, with 
an accretion efficiency rj > 0.09. 

Such a tight upper limit on the radiative efficiency can 
be somewhat relaxed if: (a) the local SMBH mass density 
has been overestimated (see the discussion in Merloni et al. 
2004); (b) the Compton Thick population has a very differ- 
ent redshift distribution from that assumed in the Gilli et al. 
(2007) XRB synthesis models (with more heavily obscured 
high luminosity objects at high redshift); (c) the primordial BH 
seeds are indeed very massive, resulting from the direct col- 
lapse of supermassive stars (Koushiappas, Bullock &. Deke l 2004[ 
[Begelman, Volonteri & Rees 2006 , Lodato & Natarayan 2007^ 

On the other hand, the strict lower limit on the efficiency could 
also be relaxed if the local SMBH mass density had been underes- 
timated, or if there was a substantial population of "wandering" 
SMBH, i.e. objects that have been kicked out of their galactic host, 
most likely for the recoil experienced due to asymmetric gravi- 
tational wave emission subsequent to a binary black hole merger 
event, after having accreted substantial amount of matter. Those 
black holes would then be computed in the AGN luminosity func- 
tions, but not among the local relic population that contributes to 

Similar calculations were shown already in Wang et al. (2008), where 
A was interpreted as the product of an "active" BH duty cycle times its 
average Eddington rate (m); with our method, that calculates the accretion 
rate distribution directly, there is no need to introduce a fixed (m) for active 
black holes. As we will discuss in ft \5A\ the very concept of active SMBH 
becomes purely observation-dependent. For a discussion of how to interpret 
lifetimes and duty-cycles see also M04. 
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Figure 4. Redshift evolution of the accretion rate density onto supermassive 
black holes 'I'bhahX^)- Darker shaded area (solid lines) shows the total 
accretion rate, while the lighter ones (dashed lines) shows the mass accreted 
in the LK mode only. This corresponds to about 20-27% of the total mass 
growth of the SMBH population. 



Pbh,o- a discussion of the relevance of this effect on the mass func- 
tion of SMBH can be found in Volonteri (2007). 

Summarizing, we can define: (i) a normalized local SMBH 
mass density = Pbh,o/(4.3 x IO^MqMpc"^); (ii) the nor- 
malized SMBH mass density at z = Zi (where Zi is our current 
observational horizon, i.e. the maximum redshift at which AGN lu- 
minosity functions are well known) = pbh{z = Zi)/PBH,o and 
(iii) the normalized mass density of SMBH ejected from galactic 
nuclei due to the gravitational wave recoil after mergers, (^lost = 
pBH,lost /pBH,o- Then, our calculations show that the average ra- 
diative efficiency is constrained to be 



0.065 



Co(l+Clost) 



< 



^rad , 



0.070 



Co(l-?i+Clost) 



(14) 



Finally, we can calculate explicitly how much mass have 
SMBH accumulated in different modes. Specifically, we find that 
between about 73% and 80% (depending on the choice of the lumi- 
nosity functions) of the relic black hole mass density has been ac- 
cumulated in radiatively efficient modes (either HR or HK), while 
only 20% to 27% in the low radiatively efficient (LK) mode, even 
if most of the time in the life of a SMBH is likely spent in this 
latter mode (Cao 2007 ; Hopki ns^ Narayan & Hemquist 2006[ >. The 
accretion rate density evolution for both the entire population and 
for the SMBH in the LK mode is shown in Figure|4] 



5.2 The SMBH mass function evolution 

So far, we have only discussed integral constraints, i.e. those ob- 
tained from the study of the evolution of the integral of the mass 
(or luminosity) function. However, the solution of eq. ^ provides 
us with the full evolution of the SMBH mass function, and we pro- 
ceed now to discuss it in detail. 

The SMBH mass function evolution is shown in Figure[5]for 



Figure 6. Differential distribution of the Eddington ratio A = ibol/^Edd 
at different redshifts. In all panels, the distribution is shown as shaded areas, 
while dashed lines represent the z = 0.1 case, shown as a reference. 



our fiducial case of erad = 0.07. In the different panels, each cor- 
responding to a different redshift, the dashed lines show the low 
redshift {z = 0.1) mass function as a reference, while the shaded 
areas represent the calculated mass functions with their own uncer- 
tainties. 

The so-called "downsizing" is apparent in the first four or five 
panels (i.e. up to redshift ~ 1.5 — 2): we clearly see how in this 
time interval the mass function grows much more substantially at 
the low-, rather than at the high-mass end. This is to be expected, in 
the light of the results of Heckman et al. (2004). There, on the basis 
of a sample of about 23,000 SDSS type II AGN for which both an 
estimate of the black hole mass (via a measurement of the velocity 
dispersion of the host) and of the intrinsic accretion rate (via mea- 
sures of the narrow line luminosity) were available, it was shown 
that it is only the small mass black holes that are growing rapidly 
at 2 ~ 0.2. Similar trends, indicating higher Eddington ratios (and 
smaller growth times, see below) have been also found recently in 
the SDSS quasar sample ( Netzer & Trakhtenbrot 2 007 ), dominated 
by un-obscured, luminous AGN at higher redshift. We stress, how- 
ever, that our results, based on a much more statistically complete 
census of the AGN population (thanks to the hard X-ray selection), 
is much less prone than previous ones to selection biases. 

Differently from what previously found in M04, we note here 
that by the epoch of the peak of the AGN/QSO activity (z ~ 2), 
the high-end of the mass function decreases significantly, too, in 
coincidence with the observed drop of the number density of bright 
QSOs at those redshifts ( IFan et al. 20011 Wall et al. 2003J 



Hasinger, Miyaji & Schmidt 2005 
Silverman et al. 2008IIBrusa et al. 2008^ . 



[Richards et al. 20061 



5.3 The accretion rate distribution function: downsizing and 
its reversal 

As we have stressed a few times already, our calculations do not as- 
sume any particular distribution of accretion rates for AGN, rather 
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Figure 5. Redshift evolution of the SMBH mass function. Sliaded areas represent the uncertainties calculated based on the estimated uncertainties in the 
luminosity functions entering our calculations (see section [5] for details). The continuity equation for SMBH growth was integrated assuming a constant 
accretion efficiency of = 0.083, which for the adopted model of the various accretion modes (§|2) corresponds to a average radiative efficiency of ~ 0.07. 
In all panels the z = 0.1 mass function is shown as a dashed line for reference. Black hole masses are measured in units of solar masses. 



derive it consistently from the observed mass and luminosity func- 
tions (see § [Sj. It is thus interesting to show the evolution of the 
(Eddington-scaled) accretion rate function, here defined as: 

$A= J <i>j,,^j^i[M,M{M,X)]^dlogM, (15) 

where (A/, M) is given by eq. l|5j, the relationship between 

M and A is simply given by A = eradAfc^/I/Edd, and A is related 
to the accretion rate via eq. (TJ. We show the Eddington ratio dis- 
tribution function evolution in Figure |6l The various panels make 
evident the change in the shape of the distribution, which results 
from the combined evolution of the X-ray luminosity function (a 
form of luminosity dependent density evolution is noticeable in the 
evolution of $a as well) and of the mass function itself. Overall, 
the trend is for a progressive flattening of the accretion rate dis- 
tribution function with increasing redshift, i.e. of a an increasing 



relative importance of highly accreting objects, as it is reasonable 
to expect on physical grounds if the SMBH population grows from 
small mass high-redshift seeds. 

The evolution of the Eddington ratio distribution function 
shown in Figure [6] does not contain any information about the typ- 
ical accretion rate as a function of black hole mass. A useful way 
to present this information is that of showing the specific growth 
rate as a function of SMBH mass. We have computed this quantity 
by taking the ratio of the black hole mass to the average accretion 
rate {M{M)). Such a ratio M/{M) defines a timescale, the so- 
called growth time, or mass doubling time, as it measures the time 
it would take a black hole of mass M to double its mass if accreting 
at the currently measured rate of M. The redshift evolution of the 
growth time (in years) as a function of black hole mass is shown in 
Figure|7] As a reference, we show in each panel the age of the Uni- 
verse as a horizontal dashed line. Black holes with growth times 
longer than the age of the Universe are not experiencing a major 
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Figure 7. Redshift evolution of the growth time (in years) as a function of black hole mass. In each panel, a horizontal dashed line marks the age of the 
Universe at that particular redshift. Black hole masses are measured in units of solar masses. 



growth phase, which must have necessarily happened in the past. 
On the contrary, objects with growth times shorter than that are ac- 
tively growing. 

The 2 = 0.1 growth time distribution is in good agreement, 
both in shape and normalization, with that first measured in the 
above mentioned work of Heckman et al. (2004). The positive 
slope of the average growth time vs. M is one of the most di- 
rect evidence of a SMBH "downsizing", i.e. of the fact that in 
the local universe only small mass black holes are actively grow- 
ing. Interestingly enough, a similar trend has been since unveiled 
for star forming galaxies jFeulner et al. 2005] INoeske et al. 20071 
IPerez-Gonzalez et al. 2008] l. where the ratio of galaxy mass Mgai 
to star formation rate is studied as a function of Mgai. It is clear that 
a direct comparison of the SMBH growth times so measured with 
those of the stars in galaxies of various morphologies/ages/masses 
may hold the key for a thorough understanding of the SMBH- 
galaxies co-evolution, and we aim at presenting such a comparison 
in a following paper. 



The redshift evolution of the growth times distribution can be 
used to identify the epochs when black holes of different masses 
grow the larger fraction of their mass. Figure [7] shows that as 
we approach the peak of the black hole accretion rate density 
{z ~ 1.5 — 2, see fig. |2j, we witness the rapid growth of high 
mass black holes, too. The growth time distribution as a function of 
mass becomes flat there, and sinks below the corresponding age of 
the Universe. This is the period of rapid, widespread growth of the 
entire SMBH population. Correspondingly, above this redshift, the 
"downsizing" trend seems to disappear. 

Finally, we show in Figure [8] the redshift evolution of the 
radiated (bolometric) energy density split into various black hole 
mass bins. On the left, we divide objects according to their mass 
at each redshift. The low redshift downsizing trend is here appar- 
ent in the rapid decline of the radiated energy density of the most 
massive black holes (log M > 9), which becomes less dramatic 
when less massive holes are considered. While in the local Uni- 
verse (i.e. for z < 0.3) the larger contributors to the total radiative 



A synthesis model for AGN evolution 1 3 



energy density are the numerous, small, rapidly accreting SMBH 
(with log M < 7), above redshift one most of the radiative en- 
ergy is produced by larger black holes, typically with masses be- 
tween 10* and 10' solar masses. As we move to higher redshift, 
the numbers of black holes with log AI > 8 declines, and the ra- 
diative energy density is dominated again by smaller black holes 
with 8 > log M > 7. The most massive black holes, those with 
log M > 9, never dominate the radiative energy density injection, 
essentially because of their small numbers, being always in the ex- 
ponentially declining part of the mass function, well above its knee 
(see Figure|5}. 

As we discussed in the Introduction, cosmological studies 
of SMBH evolution offer the advantage, with respect to those of 
galaxies, that SMBH do not undergo any morphological transfor- 
mation as they grow. Thus, once enough information is gathered 
to solve the continuity equation l[2]l, we can trace back the evo- 
lution of any black hole (or of any population of black holes of 
any given mass at a particular redshift). Indeed, it was by compar- 
ing the growth histories of black holes of different masses today 
that Marconi et al. (2004) first noticed what they called the "anti- 
hierarchical" evolution of the SMBH population (see also M04; 
Shankar et al. 2008b). Obviously, we can do the same here. On 
the right hand side of Figure [8] we show the redshift evolution of 
the radiative (bolometric) energy density split into bins of different 
black hole masses today, i.e. at 2: = 0. 

From this we see that the progenitors of local SMBH with 
masses in the range 8 < log Ai < 9 dominate the AGN radiative 
energy output between z ~ 0.3 and 2: ~ 3, when the progenitors 
of the largest black holes today (those with log M > 9) started to 
become at least equally important. 
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Figure 9. Mass functions of active SMBH with observed 2-10 keV flux 
above 10~" (red), 10~"'3 (blue) and (green) erg/s/cm^. As 

a reference, the total SMBH mass function is also shown as a yellow area 
between solid Unes, from fig.|5] 



5.4 Selection effects: mass function of "active" black holes 

From our discussion so far, it should be clear that a complete pic- 
ture of the growth and evolution of the SMBH population cannot 
be achieved without appreciating the full extent of the wide distri- 
bution of both masses and accretion rates, a point often overlooked 
in either semi-analytic or numerical models. As a consequence, the 
very concept of a clear distinction between "active" and "inactive" 
black holes can be misleading, once we consider within the same 
scheme black holes accreting close to the Eddington limit and those 
with an Eddington-scaled accretion rate as small as, say, 10^^. 

Thus, we advocate here more operational definitions of "ac- 
tive" black hole, the easiest and most straightforward of which are 
those based on nuclear flux limits. As an illustrative example, we 
plot in Figure |9] the mass functions of SMBH whose observed X- 
ray flux exceeds some fixed limits. In particular, alongside the total 
mass function, reproduced here from Fig. |5] as a term of compar- 
ison, we plot the mass functions of all SMBH with an observed 
(i.e. absorbed) 2-10 keV flux above lO"^'^ (red), lO"^*'^ (blue) 
and lo^i^ -^^ (green) erg/s/cm^. The latter two limits have been 
chosen to match the observational limits of recent deep (Chandra 
Deep Field South, CDFS; Giacconi et al. 2002) and medium-deep 
(XMM-COSMOS, XCosmos; Hasinger et al. ) surveys. Due to the 
relative steepness of the accretion rate functions at low redshifts, 
the mass functions of active black holes so defined are extremely 
sensitive to the flux limit defining what "active" means. 

This fact implies that, when trying to determine the SMBH 
mass function from a survey that selects AGN above a given flux 
limit, one should carefully account for the bias introduced by the 
fact that SMBH of any mass always have a broad distribution of 
accretion rates l IBabic et al. 2007 1 . For obvious reasons of simplic- 



ity and coherence, we discuss here hard X-ray flux limits, as we 
use the hard X-ray luminosity function to unveil the SMBH evo- 
lution, and the X-ray background to constrain the relative numbers 
of obscured and unobscured sources. However, it is clear that, pro- 
vided robust bolometric corrections are available, similar calcula- 
tions can be converted into flux limits in any band. For example, 
Greene and Ho (2007) have recently presented the mass function 
of active black holes, identified as broad Ha line emitters (in the 
redshift range 0.1 <2<0.3) in the SDSS database. Comparison 
with the soft X-ray selected AGN luminosity function of Hasinger 
et al. (2005) showed that their sample approach completeness at 
luminosity corresponding to a (rough) X-ray flux limit of 10"^^ 
erg/s/cm^. Indeed, the mass function we show at low z in Fig.|9]for 
such a flux limit is in quite good agreement with that reported in 
Greene and Ho (2007), once the fact that our active mass functions 
do include both obscured and unobscured objects is accounted for. 

In Figure [To] we show the ratio of the various "active" mass 
functions to the total SMBH mass function. These curves can be 
simply interpreted as "selection functions" for the black hole mass 
of AGN surveys of corresponding hard X-rays flux limits. However, 
as far as we can claim that, at all redshifts considered here, there is a 
one-to-one correspondence between SMBH and galaxies mass (as 
it was implicitly assumed to derive the z = mass function), then 
the various curves of Figure [TO] could be also interpreted as "AGN 
fractions" as a function of SMBH (and galaxy) mass. Once again, 
due to the main tenet of our method of calculation, namely the fact 
that growing supermassive black holes display a very broad accre- 
tion rate distribution at any given mass, the very concept of "AGN 
fraction" is strongly dependent on the sensitivity limit that defines 
what an AGN is. In general, the AGN fraction tends to be highest at 
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Figure 8. Redshift evolution of the radiative (bolometiic) energy density of AGN. Gray shaded area (between solid lines) is the total integrated energy density, 
while orange (triple-dot-dashed), green (dot-dashed), cyan (dashed) and purple (dotted) areas are the radiative energy density emitted by black holes with 
log M > 9, 9 > log M > 8, 8 > log M > 7 and 7 > log M > 6, respectively (all in units of solar masses). The right hand plot shows the same quantities 
for objects divided according to their instantaneous mass, while the right hand one is for objected partitioned according to their final (z = 0) mass. 



the high mass end, and declines at low masses. A deep X-ray sui^ey 
like the CDFS has a quite flat 'mass sensitivity' (or mass fraction) 
at 2 ~ 1, while the shallower XMM-COSMOS at the same redshift 
starts declining already for log M < 8. At higher redshift even the 
deepest surveys are highly incomplete (in a SMBH mass sense), 
and the incompleteness rises steeply with decreasing mass. In the 
case of spectroscopic studies (see e.g. Decarli et al. 2007), where 
the definition of an active galactic nucleus relies on the capability 
of spatially and spectrally resolve AGN-excited emission lines, the 
AGN fraction will be a more complicated function of SMBH (and 
galaxy) mass than shown in Figure [Tol 



6 THE EVOLVING KINETIC LUMINOSITY FUNCTION 
OF AGN 

So far, we have focused our attention on the growth of the SMBH 
population under the assumption that, whatever the mode of ac- 
cretion of a black hole, it s possible to establish a unique corre- 
spondence between radiated (bolometric) luminosity and rate of 
change of black hole mass (see section|2ll. A choice of local black 
hole mass function, an expression for the bolometric correction, 
and a choice of a value for the accretion efficiency was all that was 
needed to extract the useful information on the SMBH growth from 
the observed X-ray luminosity function. 

However, the growth of supermassive black holes through 
mass accretion is accompanied by the release of enormous amounts 
of energy which, if it is not adverted directly into the hole (see e.g. 
Narayan & Yi 1995), can not only be radiated away, but also dis- 
posed of in kinetic form through powerful, collimated outflows or 
jets, as observed, for example, in Radio Galaxies. To reiterate the 
point made before (§|2j, it is only by assessing the relative impor- 
tance of the radiative and kinetic energy output in either "low ki- 
netic" (LK), "high radiative" (HR) and "high kinetic" (HK) modes 
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Figure 10. Ratio of SMBH mass function above a given flux limit to the 
total mass function. The various shaded areas coirespond to different ob- 
served 2-10 keV flux limits, the same as those used in figure|9] i.e. 10^^^ 
(red), IQ-i"'-^ (blue) and IQ-iS-^S (green) erg/s/cm^. 
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that a full picture of the effects of SMBH growth on the environ- 
ment will possibly emerge. 

This is important, as we expect that radiative and kinetic feed- 
back differ not only on physical grounds, in terms of coupling ef- 
ficiency with the ambient gas, but also on evolutionary grounds. 
The anti-hierarchical growth of the SMBH population discussed in 
sections [5.215.41 necessarily implies a delay between the epoch of 
Quasar (HR) dominance and that of more sedate LK mode, such 
that kinetic energy feedback plays an increasingly important role 
in the epoch of cluster formation and virialization. 

In this section, we present a method to complement the in- 
formation on the cosmological growth of the SMBH population 
described in the previous sections with an analysis of the evolving 
kinetic luminosity function of AGN. As in Merloni (2004), we will 
make use of the radio luminosity function of AGN as the main ob- 
servable for our study. However, differently from the choice made 
there, we will focus our attention on the intrinsic radio core lumi- 
nosity function. The main reason for this choice is twofold. 

On the one hand, the observed relation between BH mass, ra- 
dio and X-ray luminosity (the fundamental plane of active black 
holes, MHD03) that defines the physical state of LK mode objects 
is based on the observed (5 GHz) radio core emission, not on the 
extended one. Also, Merloni and Heinz (2007) have shown that the 
intrinsic radio core luminosity of an AGN jet is a good indicator of 
its total kinetic power, once the effects of relativistic beaming are 
taken into account. On the other hand, the physical model for the ra- 
diative emission in the core of a relativistic jet is relatively well es- 
tablished piandford & Koenigl 197"9i [Heinz & Sunyaev 2003) 1, as 
opposed to the complicated physics of radio emission and particle 
acceleration in the large scale lobes of a radio jet, where environ- 
mental effects do necessarily play an important role. Obviously, the 
price to pay for the use of radio core emission as a tracer of the jet 
power, is that relativistic beaming may severely affect the derived 
parameters, and must be taken into account. 

In the classical model by Blandford & Koenigl (1979), the 
flat spectrum radio synchrotron emission of a compact jet core is 
produced by a superposition of self-absorbed synchrotron spectra, 
each from a different region in the jet. The model predicts a de- 
pendence of the monochromatic radio luminosity on jet kinetic 

17 /12 

power VKkin of the form oc W^.^^ . More generally, Heinz & 
Sunyaev (2003) showed that any scale invariant jet model produc- 
ing a power-law synchrotron spectrum with index a„ must obey 
the relation L„ oc y/^'^'^^°'"'>/^'^ . By studying a sample of 
nearby low-luminosity radio AGN for which the kinetic power had 
been measured, Merloni & Heinz (2007) have shown that, once 
beaming is statistically taken into account, a tight relation between 
the 5 GHz intrinsic radio core luminosity (Lr) and the jet kinetic 
power holds, in the form: 



erg s 



(16) 



where Wo ~ 1.6 x 10^*^ and Lq = 10^° erg s"\ and the intrin- 
sic scatter is about 0.35 dex. In i]6.1l we will use this relation to 
construct the kinetic jet luminosity function from the observed flat 
spectrum radio luminosity function ^l{Lv) (abbreviated as FSLF 
below). Similar, but complementary studies of the AGN kinetic 
luminosity function evolution based on steep spectrum sources 
has been recently carried out jKording, Jester & Fender 2008 [ 
IShankar et al. 2008at . The results of both these works are in rea- 
sonable agreement with those we present here. 



6.1 Beaming effects and the intrinsic luminosity function of 
radio jet cores 

As pointed out in Heinz, Merloni & Schwab (2007), we can use 
eq. M6\ and the observed FSLF to derive the underlying kinetic 
luminosity function of radio jet cores: 



(^o( 



Wo ) 0.81 Wo \Wo) 



(17) 



Without loss of generality, and consistent with most observational 
studies, we will use a broken power-law to describe (?!)r(-Lr): 



0r(-£'r) 



Lb. 



Lr 



(18) 



In the above equation il6\ . however, it is the intrinsic radio core 
luminosity that enters the determination of the kinetic power Since 
jets are relativistic, the observed luminosity functions, 0' are al- 
ways affected by Doppler boosting. Assuming all object have the 
same bulk Lorentz factor 7 = (1 — the observed lumi- 

nosity L of a relativistic jet is related to the emitted one, C, via the 
relation L — S^C, where S, the kinematic Doppler factor of the 
jet is defined by 5 = [7(1 — f3cMs6)]~^. The exponent p depends 
on the assumption about the jet structure, and is here taken to be 
p = 2 as appropriate for continuous jets, rather than discrete ejec- 
tions (for which p = 3, see Urry & Padovani 1995). If the radio 
sources have all a pair of oppositely directed jets, and the jets are 
randomly oriented in the sky, the cosine of the angle between the 
velocity vector and the line of sight, cos have a uniform proba- 
bility distribution, and the observed sources will be those lying at 
angles between 6'min and Smax, corresponding to Doppler factors 
of <5min and (5max, respectively (see Urry and Shafer 1984). We as- 
sume here that ^max ~ 0° and Smin = I/7 (i.e. neglecting all 
strongly de-beamed objects). 

Urry and Padovani (1991) have shown how, under the above 
assumptions, an intrinsic double power-law LF is transformed by 
the effects of relativistic beaming. It can be shown that, for astro- 
physically relevant values of ai, 02 and 7 the most general sit- 
uation is that in which the beamed LF is a triple broken power- 
law ( [Heinz, Merloni & Schwab 2007) . The faint- and bright-end 
slopes are unchanged, while around the knee an intermediate-slope 
-{p+l)/p= -3/2 is found. 

The problem of inverting the Urry and Padovani (1991) so- 
lution, i.e. of deriving the intrinsic shape of the radio loud AGN 
LF given the observed one and a few assumption about p and 7 is 
greatly simplified if the intermediate part of the beamed luminosity 
function is neglected. In practice, in such a case the intrinsic, de- 
beamed, FSLF <j) can be obtained from the observed one simply by 
the following transformations of the LF knee, while keeping both 
end slopes fixed. In detail, following the procedure described in 
Urry and Padovani (1991), the LF normalization will be corrected 
by a factor 



-l/(a 



(19) 



while the break luminosity transforms as: 

Lk,c(^) ^ Ln,4z) = Lk,c(^)(^C71/C2)-i/(''^-'^i\ (20) 

where ki = S^J, - S^^, K2 = <5giax - Chf' = /m and 
CI = ai - - 1, C2 = aa - (1/p) - 1. 
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Figure 11. Kinetic luminosity function of AGN jets (here plotted as W\^i^ X (jtw) at four different redshifts. In each panel the green area between the black 
solid lines represents the total kinetic luminosity function, with the uncertainty reflecting the difference in the luminosity functions between Dunlop & Peacock 
(1990) and De Zotti et al. (2005). The blue area between blue dashed lines represents the contribution from sources in the LK mode, while the dotted lines 
mark the contribution from sources in the HK state. In the Left plot (FS) the jet core luminosity function includes only flat spectrum sources, while in the right 
plot (FS+SSmax) a maximal contribution from hidden cores in steep spectrum sources is included (see text for details). 



In this way, for any given luminosity function and mean 
Lorentz factor, it is possible to derive the intrinsic kinetic lumi- 
nosity function of the AGN jet, by simply applying eq. (TT} to the 
de-beamed radio luminosity function. 

As opposed to the case of X-ray selected luminosity func- 
tions, the available constraints on the radio cores luminosity func- 
tions are not very stringent. Here we will use two different deter- 
minations, translated into our Cosmology, and computed at 5GHz: 
the flat spectrum radio luminosity function from Dunlop & Pea- 
cock (1990) (their HIGH-z parametrization of their table C4), and 
that from De Zotti et al. (2005), given by the sum of flat spec- 
trum radio quasar and BL Lac objects as in eqs. (l)-(3) and ta- 
ble 1 of their paper. As we will show in the following, and as al- 
ready discussed in Heinz, Merloni and Schwab (2007), the total ki- 
netic power released by AGN is sensitive to the faint-end slope of 
the kinetic luminosity function. For consistency with observations 
of local radio luminosity functions down to very low luminosities 
( |Nagar, Falcke & Wilson "2005] >, we have fixed the 2 = faint end 
slope to ai = 1.85 for both the FSLF of Dunlop & Peacock (1990) 
and the BL Lac population of De Zotti et al. (2005). Because ai — 1 
is (marginally) larger than 0.81, the slope of the Wkin-Ln relation 
( 116) , the total kinetic power is determined by the lower cut-off of 
our radio luminosity function, chosen here to correspond to the ra- 
dio luminosity above which the total number of sources is equal 
to the total number of SMBH with mass above 5 x IO^A^q (see 
§|4j4j. However, the difference between the total kinetic power ob- 
tained in this way and that one would get for oi < 1.81 is smaller 
than the uncertainty introduced by the intrinsic scatter in the re- 
lation ( |16t . Therefore, for the ease of computation, we have al- 
lowed a very mild evolution (flattening) of such a faint-end slope 
between 2 = and 2 = 2, where we have fixed ai = 1.75. In- 
deed, there are independent lines of evidence pointing towards a 
flattening of the faint-end of the AGN radio luminosity function at 



high redshift ( jCirasuolo, Magliocchetti & Celotti 2005) , similarly 
to the well known evolution observed in the X-ray band. Future 
better constraints on the high redshift evolution of the faint radio 
AGN population will surely help tightening our constraints on the 
overall kinetic power density. 

The observed luminosity functions are then de-beamed fol- 
lowing the above described prescription, assuming p = 2. For the 
mean jet Lorentz factor, we adopt the most likely value based on 
the statistical analysis of core radio luminosity of a small sample of 
15 AGN with measured kinetic power (see Merloni and Heinz 2007 
for detailsj^: 7 = 7. Also, for the sake of simplicity, and lacking 
any clear observational evidence, we do not account here for the 
possibility of varying Lorentz factors with redshift. 

So far, we have proceeded under the assumption that the ob- 
served FSLF are complete, in the sense that they account for the 
totality of jet cores. In fact, we have to discuss also the possible 
contribution from flat spectrum cores which are 'hidden' below 
a more powerful, extended steep-spectrum component. By defini- 
tion, the FSLF contains all black holes except those included in the 
steep-spectrum luminosity function (SSLF), which are dominated 
by optically thin synchrotron emission. The scaling relation from 
eq. l ll6t does not hold for extended, steep-spectrum sources (but 
see Birzan et al. 2004; Best et al. 2006). We can, however, derive 



* We note here that the majority of studies of jet core velocity struc- 
tures at VLBI scales have recently shown that the distribution of intrin- 
sic Lorentz factors is relatively broad, typically in the range 3 < 7 < 30, 
jLaing et al. 1999|lCohen et al. 20071 . Although cumbersome, it is straight- 
forward to compute the de-beamed radio luminosity function assuming 
more complicated distributions of the bulk Lorentz factors. However, for 
the sake of clarity in the current discussion, we leave this exercise to further 
studies. 
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an upper limit to the contribution from steep spectrum sources. As- 
suming a typical optically thin synchrotron spectral index of 0.65, 
any underlying flat spectrum component would have to fall be- 
low ~20% of the observed 5GHz luminosity, otherwise the source 
would become too flat to qualify as a steep-spectrum source. Us- 
ing the steep-spectrum luminosity functions (SSLF) from Dunlop 
& Peacock (1990; high-z parametrization of their table C4) and De 
Zotti et al. (2005), we then calculate a maximal core luminosity 
function comprising both flat spectrum sources and those 'hidden' 
cores for which we take a contribution equal to 10% of the SSLF. 
This can be considered as an extreme case in which the objects 
with powerful (optically thin) radio lobes have a core-to-lobe flux 
ratio (the so-called core prominence, Bridle et al. 1994) of 0.1 (for 
a discussion on the observational constraints on the intrinsic core 
prominence, see Laing et al. 1999, Jackson and Wall 1999). 

The kinetic luminosity functions calculated from the observed 
radio cores luminosity function with the help of eq. l |16b are shown 
in Figure [TT] where we plot the luminosity function multi- 
plied by the kinetic power Wkin in a representation that highlights 
the luminosity range where most of the power is released. For the 
sake of clarity we have chosen to keep separate the results obtained 
from adopting only the FSLF (left panel [FS]) from those obtained 
adding to the FSLF a contribution of 10% of the observed SSLF 
(right plot [FS-l-SSmax]). As discussed above, we regard this latter 
as an effective upper limit to the AGN kinetic luminosity function. 
Figure nn also shows the decomposition of the total kinetic lumi- 
nosity functions into the two kinetically dominated modes of ac- 
cretion, as discussed in § [2] For both cases (FS and FS-l-SSmax), 
at z ^ 1, the kinetic luminosity functions is dominated by low lu- 
minosity (LK) objects, with bright, radio loud QSOs (HK) only 
contributing to the high power tail. At, and above, the peak of the 
BHAR density evolution {z ^ 1.5 — 2.5, see Fig. |2l HK objects 
start dominating the total kinetic power output of the growing black 
holes population, even if the relative contribution of LK sources is 
always substantial. 



6.2 Kinetic energy density evolution and kinetic efficiency 

Integrating the kinetic luminosity functions obtained following the 
procedure described in the previous section, we arrive to an esti- 
mate of the total kinetic energy density as a function of redshift. 

Figure [T2] shows the redshift evolution of the total integrated 
kinetic energy density (divided by ric?, so that it can be displayed in 
units of Mq yr~^ Mpc^'^). Once more, we show separately the re- 
sults of the calculation done including flat spectrum sources only 
(left panel) and one in which a maximal contribution from hid- 
den cores in steep spectrum sources is added (right panel). The to- 
tal jet kinetic power density (shown as grey shaded areas) is split 
into the contributions from LK and HK modes of growth. Because 
the average Eddington ratio increases towards high redshift, so it 
does the relative contribution of HK objects (i.e. radio loud QSOs, 
or powerful FR II), irrespective of the choice of core radio lumi- 
nosity function. As a term of comparison, the plot shows also the 
evolution of the total black hole accretion rate density (BHAR), 
which we discussed at length in section 15.11 Also shown is the 
estimated total kinetic power output of SNe II, as computed di- 
rectly from the best fit to the star formation rate density evolution 
( [Hopkins and Beacom 2006) . The total kinetic power released by 
low luminosity (LK) AGN at z < 0.3, where it is by far the domi- 
nant contributor to the total jet kinetic power, is comparable to that 
released by type II supemovae, which, however, become about a 



factor of few (up to about an order of magnitude) higher at the peak 
of star formation activity. 

Given the stark differences in the redshift evolution of low- 
and high-power sources, and the fact that the relative importance 
of kinetic power output increases as the Eddington ratio decreases 
(see §|2]l, we can conclude that the most important feature of the 
evolution of the kinetic energy output of growing black holes, as 
compared to the radiative power output (see Figure [8]( is the very 
mild redshift evolution. In this respect, it is interesting to mention 
that Cavaliere & Lapi (2008) have recently shown how the problem 
of the 'missing baryons' in groups and cluster of galaxies can be 
understood in terms of a 'dual' model for AGN feedback (kinetic 
at low z and radiative at high z) that is in very good agreement 
with the overall picture of SMBH growth and energy output that 
emerges from the work presented here. 

A fundamental question we are seeking the answer for is the 
accurate determination of the kinetic energy production efficiency 
of growing black holes. This bears obvious important consequences 
for AGN feedback models, but has so far been treated as a free 
parameter. Our calculations allow for a more precise determination 
of such parameter. Let us first compute the total integrated (mass 
weighted) average kinetic efficiency as: 

(6k- > = /o"' Pkin(^)rfz 

^ f 0.0028 ± 0.0008 FS 

~ \ 0.0053 ± 0.0015 FS + SSmax 

These numbers are consistent with those found in Heinz, Merloni 
and Schwab (2007), once we consider the mean Lorentz factor of 
7 = 7 used to calculate the intrinsic FSLF, and the steeper slope of 
the I/R-Wkin correlation l ll6t . On the other hand, they are some- 
what smaller (less than a factor of 2) than those found in both 
Kording, Jester & Fender (2008) and in Shankar et al. (2008a), 
where, however, only steep-spectrum sources were used, and the 
total kinetic power was derived from the correlation between ex- 
tended radio luminosity and kinetic power of Willott et al. (1999). 

Combining the results discussed in section [STI on the aver- 
age radiative efficiency with those shown in eq. | |2U we conclude 
that SMBH during their growth from z ~ 5 till now convert about 
15-25 times more rest mass energy into radiation than into kinetic 
power, with the exact number depending on the poorly known de- 
tails of the intrinsic jet cores luminosity function, as well as on 
our assumptions about the beaming corrections to be made. We 
note here that numerical simulations of AGN feedback in merging 
galaxies ( |Di Matteo, Springel & Hernquist 2005| l fix the feedback 
efficiency to about 5 x 10"'^ in order to reproduce the M — cr re- 
lation. In light of what we have shown, it could be argued that this 
feedback might be almost entirely provided by kinetically domi- 
nated modes of SMBH growth, as long as the coupling efficiency 
between the kinetic power provided by AGN and the interstellar 
and intra-cluster medium is equal to one. 

We can also compute directly the kinetic efficiency as a func- 
tion of redshift and SMBH mass, which we show in Fig.[T3] In each 
of the two panels, corresponding to the FS and FS-FSSmax cases, 
the horizontal solid lines show the mass weighted average from 
eq. The various curves describe the main properties of kinetic 
feedback as we observe it. For each of the chosen mass ranges, the 
kinetic efficiency has a minimum when black holes of that mass 
experience their fastest growth: this is a different way of restating 
the conclusion that most of the growth of a SMBH happens during 
radiatively efficient phases of accretion. However, when the mass 
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Figure 12. Redshift evolution of the kinetic energy density (divided by rjc^ and plotted in units of Mq yr~^ Mpc^"') is shown with dotted lines and gray 
area. The top solid lines show the BHAR density evolution. 'I'bharCz) (see fig.|4). while red (dashed lines) and blue (dot-dashed lines) bands represent the 
relative contribution to the total kinetic power from LK and HK sources, respectively. The green solid line shows the estimated total kinetic power output from 
SNe II, as inferred by Hopkins and Beacom (2006) based on the measured star formation rate density. The left plot shows the case in which only fiat spectrum 
radio AGN are include in the intrinsic core radio luminosity function (FS), while the right hand plot shows the case in which a maximal contribution from 
'hidden' flat spectrum cores in steep-spectrum sources is included. 
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Figure 13. Redshift evolution of the kinetic efficiency e^in. SMBH of different masses are here plotted separately, with a color coding analogous to that of 
Fig lH The right hand plot shows the results obtained including flat spectrum radio sources only, while the right hand one shows the results obtained including 
also a maximal contribution from the hidden cores of steep spectrum radio AGN. In each plot, horizontal black solid lines mark the mass weighted average 
values for the kinetic efficiencies, with the dashed lines representing the uncertainties from the particular choice of radio and X-ray luminosity functions. 
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increases, SMBH are more and more likely to enter the LK mode 
(see section[2ll, which increases their kinetic efficiency. More mas- 
sive holes enter this phase earlier, and by 2: = they have reached 
the highest kinetic efficiency of ~ 2 5 x 10~^, as it is clearly 
seen in Figure [T3] This is a natural consequence of the observed 
anti-hierarchical growth of the SMBH population, and of the cho- 
sen physical model for the accretion mode of low-Eddington ratio 
objects. 



7 DISCUSSION AND CONCLUSIONS 

We have presented a comprehensive review of our current knowl- 
edge of the evolution of AGN and of the associated growth of the 
supermassive black holes population. 

Similar to the case of X-ray background synthesis models, 
where accurate determinations of the XRB intensity and spectral 
shape, coupled with the resolution of this radiation into individ- 
ual sources, allow very sensitive tests of how the AGN luminosity 
and obscuration evolve with redshift, we have argued that accurate 
determinations of the local SMBH mass density and of the AGN 
(bolometric) luminosity functions, coupled with accretion models 
that specify how the observed AGN radiation translates into a black 
hole growth rate, allow sensitive tests of how the SMBH popula- 
tion (its mass function) evolves with redshift. By analogy, we have 
named this exercises 'AGN synthesis modelling' . In performing it, 
we have taken advantage of the fact that the cosmological evolution 
of SMBH is markedly simpler than that of their host galaxies, as 
individual black hole masses can only grow with time, and SMBH 
do not transform into something else as they grow. Moreover, by 
identifying active AGN phases with phases of black holes growth, 
we can follow the evolution of the population by solving a simple 
continuity equation (|2), where the mass function of SMBH at any 
given time can be used to predict that at any other time, provided 
the distribution of accretion rates as a function of black hole mass 
is known (see section[3]l. 

Such a straightforward approach is fully justified as long as 
SMBH coalescences, which are likely to occur in the nuclei of 
merged galaxies, do not alter significantly the shape of the mass 
functions. It is well established, as we have also shown here, that 
the local black hole mass density is fully consistent with the mass 
accumulated onto SMBH in AGN episodes for reasonable values 
of the accretion efficiency (see § 15. It . Mergers only redistribute 
mass so they do not affect any integral constrain, such as the Soltan 
argument, and variations thereof. Until now, evaluating the rela- 
tive importance of SMBH mergers in the mass function evolution 
has proved itself very difficult dYu & Tremaine 2002t . due to the 
large uncertainties in both merging rates of galaxies (either theoret- 
ically or observationally, see e.g. Bundy, Treu & Ellis 2007; Guo 
and White 2008; Hopkins et al. 2008 and references therein) and 
in the physical processes involved in the merger of the two nuclear 
black holes (see Colpi et al. 2008 for a recent review). Shankar et 
al. (2008b) have attempted to estimate such effect considering only 
equal mass merger, at a rate most likely far in excess of the true 
one, and concluded that the effect of SMBH mergers on the local 
mass function is overall smaller than the current uncertainties in the 
AGN bolometric luminosity function itself, and may be relevant, if 
at all, only at the very high mass end of the distribution. A more 
detailed comparison, taking into account the full mass and redshift 
dependencies of the merger rates, is however necessary before a 
firm conclusion can be drawn. 

In order to carry out our calculation, we assumed that black 



holes accrete in just three distinct physical states, or "modes": a ra- 
diatively inefficient, kinetically dominated mode at low Eddington 
ratios (LK, the so-called "radio mode" of the recent literature), and 
two modes at high Eddington ratios: a purely radiative one (radio 
quiet, HR), and a kinetic (radio loud, HK) mode, with the former 
outnumbering the latter by about a factor of 10. Such a classifi- 
cation is based on our current knowledge of state transitions in 
stellar mass black hole X-ray binaries as well as on a substantial 
body of works on scaling relations in nearby SMBH (see § |2j. It 
allows a relatively simple mapping of the observed luminosities 
(radio cores. X-ray and/or bolometric) into the physical quantities 
related to any growing black hole: its accretion rate and the released 
kinetic power. 

With the aid of such a classification, we have solved the con- 
tinuity equation for the black hole mass function using the lo- 
cally determined one as a boundary condition, and the hard X- 
ray luminosity function as tracer of AGN growth rate distribu- 
tion, supplemented with luminosity-dependent bolometric correc- 
tions tMarconi et al. 2004) and absorbing column density distribu- 
tions. By construction, our adopted intrinsic hard X-ray luminosity 
functions satisfy the XRB constrain (as well as the sources number 
counts, and many others), as they adopt the most recent synthesis 
model as a input (see Gilli et al. 2007 for details). One remain- 
ing uncertainty with these models lies in the actual luminosity and 
redshift distribution of heavily obscured (Compton thick) AGN. 
However, according to the adopted column density distribution, the 
contribution to the overall mass growth from Compton thick AGN 
amounts to just about 20%. Firmer conclusions on this issue will 
need more accurate determination of number density and typical 
luminosity of these heavily absorbed (and thus invisible in the 2-10 
keV band) sources at redshift between 1 and 2. In fact, preliminary 
results based on mid-infrared selection criteria (Daddi et al. 2007] 
IFiore et al. 2008l l seem to show quite a good agreement with the 
predictions of the XRB synthesis models adopted here. 

The main results of our study are summarized below. 

• We have shown how using standard Soltan (1982) type of 
arguments, i.e. comparing the local mass density to the inte- 
grated mass growth in AGN phases, very tight constraints can be 
put on the average radiative efficiency of the accretion process: 
eo(i+Lo ~ < where is the local mass den- 

sity in units of 4.3 xIO^Mq Mpc~^, while and ^lost are the 
mass density of 2; 5 and "wandering" SMBH, respectively (also 
in units of the local mass density). Figure[74l shows the constraints 
derived from eq. l lI4b in the ^lost, (a*) panel, where (a*) is the 
mass weighted average spin parameter of the SMBH calculated in- 
verting the classical GR 77(a) relation (Shapiro 2005) and taking 
into account the assumed relationship of eq. (T} between iq and 
trad- It is interesting to note that both the amount (numbers and 
masses) of black holes effectively ejected from galactic nuclei due 
to gravitational wave recoil after a merger and the average radiative 
efficiency of accreting black holes depend on the spin distribution 
of evolving SMBH (see e.g. Volonteri 2007, Berti and Volonteri 
2008). Thus, eq. ( I14t couples implicitly the spin distribution of ac- 
creting black holes (through ^lost and trad) and the properties of 
the seed black hole population, whose density must be reflected in 
the 2 ~ 5 mass density ^i. At face value, our results indicate that, 
for ^0 = 1 3> Ci, SMBH must have on average a rather low spin, 
similar to what predicted by models of AGN fuelling that take into 
account the self-gravity of an accretion disc and result in accretion 
events which are small and randomly-oriented, thus favoring effec- 
tive spin-down of the hole ( [King, Pringle & Hofmann 2008[ l. 
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(dot-dashed) line represent the allowed area for = 0.1 (gi = 0.3). 

• We confirm previous results and clearly demonstrate that, at 
least for z < 1.5, SMBH mass function evolves anti-hierarchically, 
i.e. the most massive holes grew earlier and faster than less mas- 
sive ones (§ 15. 3t . By looking at the distribution of SMBH growth 
times as a function of mass and redshift, for the first time we find 
hints of a reversal of such a downsizing behaviour at redshifts 
close to and above the peak of the black hole accretion rate den- 
sity (z ~ 2), where we witness the epoch of rapid growth of al- 
most the entire SMBH population. Our results bring the study of 
AGN 'downsizing' from a phenomenological level, in which one 
simply describes the observed behaviour of the luminosity func- 
tions ( |Hasinger, Miyaji & Schmidt 2005] [Bongiomo et al. 2007| >, 
to a more physical one, in which such an observed behaviour is 
explained in terms of evolution of physical quantities (mass and 
accretion rate), which could be more directly related to the parallel 
evolution of the host galaxies population. 

• Differently from most previous analytic, semi-analytic and nu- 
merical models of black holes growth, we do not assume any spe- 
cific distribution of Eddington ratios, rather we determine it empir- 
ically by coupling the mass and luminosity functions and the set of 
fundamental relations between observables in the three accretion 
modes. SMBH always show a very broad accretion rate distribu- 
tion, and we have highlighted the profound consequences of this 
fact for our understanding of observed AGN fractions in galaxies, 
as well as on the empirical determination of SMBH mass functions 
from large surveys (ij |5.4t . 

• We have presented the cosmological evolution of the kinetic 
luminosity function of AGN, based on the evolution of their flat 
spectrum radio luminosity function and on the empirical correlation 
between intrinsic radio core luminosity and kinetic power found in 
Merloni and Heinz (2007). As opposed to the mass growth evolu- 
tion, the kinetic luminosity function so derived is not very tightly 
constrained due to poor observational information on the true (in- 



trinsic) radio core luminosity function at high redshift and to the 
uncertain distribution of the AGN jet bulk Lorentz factors needed 
to correct for relativistic beaming of their core emission. Neverthe- 
less, we were able to measure the overall efficiency of SMBH in 
converting accreted rest mass energy into kinetic power, the "ki- 
netic efficiency" tkin, which ranges between 3 and 5 xl0~^, de- 
pending on the choice of the radio core luminosity function. This 
value is somewhat smaller (but less than a factor of 2) than those 
found independently in both Kording, Jester & Fender (2008) and 
in Shankar et al. (2008a), where only steep-spectrum sources were 
used, and the total kinetic power was derived from the correlation 
between extended radio luminosity and kinetic power of Willott 
et al. (1999). Nevertheless, the redshift distribution of the kinetic 
power density released by growing SMBH is in good agreement 
with that derived by Kording et al. (2008), lending further support 
to the idea that current measures of the kinetic AGN power are 
reaching quite robust conclusions. 

• As for the local (z = 0) AGN kinetic power density, pkin, we 
found it ranges between 2.5 and 10 x lO"^'"* erg s~^ Mpc"'^, compa- 
rable with the total kinetic power density from type II Supernovae 
(psNii ~ 4 X 10^^ erg s"^ Mpc~^, Hopkins and Beacom 2006) 
while the local stellar luminosity density is about p, ~ 2 x 10*^ 
erg s^^ Mpc^^ and the total AGN radiative density is about prad — 
1.6 X 10^" erg s^^ Mpc^^ 

• Our capability to resolve the mass and accretion rate functions 
allows us to separate the evolution of both growth rate and kinetic 
energy density into different mass bins and into the various modes 
of accretion. In doing so, we found that most (73-80%) of the local 
mass density was accumulated during radiatively efficient modes 
(either HR or HK), with the remaining 20-27% in the radiatively 
inefficient LK mode. The kinetic power density at low redshift is 
completely dominated by low luminosity AGN, while the contri- 
bution from radio loud QSOs (mode HK) becomes significant at 
z ^ 2 (see section \6l2\ . The measured tkin varies strongly with 
SMBH mass and redshift, being maximal for very massive holes at 
late times, a property in agreement with what required for the AGN 
feedback by many recent galaxy formation models. 

The richness of details we have been able to unveil for the cos- 
mological evolution of supermassive black holes demonstrates that 
times are ripe for comprehensive unified approaches to the multi- 
wavelength AGN phenomenology. At the same time, our results 
should serve as a stimulus for semi-analytic and numerical model- 
ers of structure formation in the Universe to consider more detailed 
physical models for the evolution of the black hole population. 
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